Ab-initio and Conformational Analysis of a Potent VEGFR-2 Inhibitor: A Case Study on Motesanib.

Vascular endothelial growth factor receptor-2 (VEGFR-2); a cell surface receptor for vascular endothelial growth factors, is a key pharmacological target involved in the cell proliferation/angiogenesis. It has been revealed that VEGFR-2 induces proliferation through activation of the extracellular signal-regulated kinases pathway. In this regard, targeting the VEGFR-2 has been considered as an efficient route to develop anti-tumor agents. Motesanib is a small-molecule antagonist of VEGFR-1, 2, and 3 (IC50s; 2 nM, 3 nM, 6 nM, respectively). It is an experimental drug candidate undergoing clinical trials against some types of cancer. In the present study, Motesanib (AMG 706) was evaluated in terms of its binding energies with individual amino acids of VEGFR-2 active site (amino acid decomposition analysis). For this purpose, functional B3LYP associated with split valence basis set using polarization functions (Def2-SVP) was used. Comparative conformational analysis of the ligand in optimized and crystallographic states revealed that Motesanib does not necessarily bind to the VEGFR-2 active site in its minimum energy conformer.


Introduction
Vascular endothelial growth factors (VEGFs) bind to three structurally related receptor tyrosine kinases; VEGFR-1, VEGFR-2, and VEGFR-3. A number of coreceptors such as neuropilins that lack intrinsic catalytic activity bind to VEGF and also modulate the effect of the VEGFRs (1). VEGFRs have a high degree of homology within the kinase domain; however, their signaling properties greatly differ (1). VEGFR-2 is the major mediator of responses in endothelial cells and it is considered to be a principal signal transducer in angiogenesis (2). Vascular endothelial growth factor receptor 2 (VEGFR-2) is a cell surface receptor for VEGFs (3,4). VEGF signaling pathway has been well demonstrated to induce angiogenesis during tumor development (5, 3). Vascular endothelial growth factor (VEGF) has been found to be a major driver of tumor angiogenesis leading to efforts in development of novel therapeutics aimed at inhibiting its activity. It is generally accepted that the destructive growth of tumors and their metastases is highly depended on angiogenesis (6).
Anti-VEGF therapy has been regarded as a prominent choice for the management of several human malignancies. The vast majority of solid comprise mainly two categories; ligandbased and structure-based methods. Within this scenario, theoretical methods based on density-functional theory (DFT) are going to play an increasingly important role in many applications of computational chemistry to drug discovery (27)(28)(29). One of the mostly used DFT methods is Becke three-parameter Lee-Yang-Parr (B3LYP) hybrid density functional theory (30). This technique has been applied for conformational analysis of antiangiogenic agents such as Motesanib (31). Estimating the proportion of individual amino acid-ligand interaction energies in total binding energy would be a very useful trend in pharmacophore discernment and development (Amino acid decomposition analysis). No such report could be found for Motesanib with VEGFR-2.
In continuation to our interest in developing new cytotoxic agents (32), we aimed to determine the binding pose and binding energies of Motesanib in VEGFR-2 active site. For this purpose, DFT calculations at the level of B3LYP/Def2-SVP were used to estimate the individual amino acid-ligand interaction energies. Furthermore; conformational analysis of Motesanib was performed at the same level of theory to determine the torsional deviation from minimum energy state upon binding to the receptor.

Experimental
Computational section X-ray crystallographic structure of VEGFR-2 with its cognate ligand (AMG 706) was downloaded from Brook Haven Protein databank (PDB code: 3EFL, www.rcsb.org). All the preprocessing procedure for ligand and receptor crystallographic files were done within WHAT IF server (European Molecular Laboratory Heidelberg, Germany). All hydrogens were properly added to the receptor PDB file using What if server. All computational DFT studies were performed using functional B3LYP associated with split valence basis set using polarization functions (Def2-SVP) (33).
The evaluated amino acid residues in ab initio study (Leu889, Ala866, Lys868, Glu885, Val899, Val916, Cys919, Leu1035, Asp1046) tumors are involved with VEGF overexpression (7). In this regard, chemical agents with the inhibitory activity on angiogenesis may be focused as a workable treatment options for patients involved with solid tumors (8). Inhibition of VEGF has been reported to significantly suppress tumor angiogenesis in mouse tumor models (9). The angiogenesis process of solid tumors may be avoided due to the inhibition of the tyrosine kinase VEGFR-2 signaling pathway. In this way required blood flow for developing tumor would decrease dramatically and even the growth will stopped because of lack of nutrient and growth factors supported by freshly forming vessels (10,11). Hence, anti-tumor drug development targeting VEGFR-2 signalling pathway has been recently highlighted as an important way in the clinical trials (12).
In lead-drug development strategies, combination of experimental methods with computer aided molecular design (CAMD) techniques is essential for the development of new drugs aimed at new targets, and thus for medicinal chemistry (26). Various computational chemistry methods are in hand for running CAMD. These methods were chosen according to the information from crystallographic data (34). The proposed interaction profile could be found in Brook Haven Protein databank. All amino acids were considered in their real electrostatic state. For each residue under study, N-terminal was acetylated and C-terminal was methyl amidated to mimic the original electron density. All conformational and configurational features were the same as the X-ray structure. However the positions of hydrogen bonds are not clearly recognized in a typical X-ray crystallographic file and this restriction persuaded us to further optimize the heavy atom-hydrogen bonds using constrained optimization method (heavy atom fixing approximation).
All ligand-receptor interaction energies were estimated by B3LYP/ Def2-SVP method. The whole calculations were done with the ORCA quantum chemistry package (30). The optimized structures were visualized by Visual Molecular Dynamics program (VMD) (35).

Ab initio study of ligand-receptor interactions
To obtain a binding profile between a Motesanib and VEGFR-2 active site, relevant amino acids were chosen on the basis of information from Protein Data Bank (PDB). The representation of the Motesanib structure in the active site of VEGFR-2 was further confirmed via schematic 2D interaction profile generated by LIGPLOT Figure 1.
The directionality of hydrogen bonds including optimum distances and angles supports efficient interactions with receptor. Bearing this in mind, the optimization process was done with the same basis set to obtain the exact geometry of H-bonds. The related data are summarized in Table 1. Hydrogen bond geometries were described as H-donor-acceptor angles. It should be noted that hydrogen bond lengths were obtained considering H-acceptor distances.
Ligand-residue binding energies (ΔE b ) were calculated using the following equation: (1) E LR stands for residue-ligand interaction energy, E R and E L indicate the electronic energies for residues and ligand, respectively. Individual ligand-residue binding energies are shown in Figure 2.
Lipophilic contacts are a function of orientation, constant/induced dipole moments and distance (36). Attractive hydrophobic interaction was made between AMG 706 and Ala866 residue (-0.56 cal/mol). Regarding the obtained results, one may assume that electrostatic hydrogen bonding interactions may have significant contribution to total binding energy between AMG 709 and receptor. For Leu889, Val899, Val916 and Leu1035 the positive binding energies might be related to an inappropriate orientation of ligand in the active site of receptor (crystallographic state). It should be noticed that molecular dynamic may assist in balancing these close contacts which are responsible for repulsive interactions with receptor.
The role of cation-pi interaction as a force for molecular recognition in biological media has been revealed via studies on model systems and the analysis of biological macromolecular structures (36). In Lys868 the quaternary amine moiety might be responsible for observed cation-π interaction with central pyridine ring of AMG 709 (-3.93 kcal/mol). Binding pose of the ligand revealed that nitrogenous cation centered on the top of the π face of pyridine ring ( Figure 3). Our estimated energy for the associated cation-π interaction correlated well with previously reported data (37), but might be less than accepted values (38). This energy difference might be attributed to the existence of pyridine lone pair. In pyridine, the lone pair does not participate in aromaticity and thus electronegativity of the heteroatom wins out and  weakens the cation-π binding ability. The deprotonated Glu885 is the most significant residue for enzyme-inhibitor interactions due to the strong hydrogen bond with the amide NH in AMG 706. Typical charge-assisted hydrogen bond between Glu885 and amide NH was found to be supported by significant binding energy of -26.82 kcal/mol. Another important H-bond between Asp1046 and amide oxygen was associated with -12.96 kcal/mol interaction energy. In the case of Cys919 residue, the observed binding energy was estimated to be -9.06 kcal/mol. In fact hydrogen bonding interaction with Cys919 from H-donor group of the inhibitor is the key feature of VEGFR inhibitors (39).
In ligand-receptor interaction, stereoelectronic effects are prominent in determining complementary potential electrostatic surfaces. Ligand electronic structure may address its proper orientation in the enzyme active site and potent inhibition would be expected regarding proper fitness of the ligand and electronic surfaces of the active site.
Mulliken partial electronic charges were assigned to the constituent atoms of compound AMG 709 (Figure 4) (40). It should be noted that atoms participated in key bindings with Asp1046, Glu885 (charge-assisted interactions) and Cyc919 residues possessed relatively negative electronic charges.

Comparative conformational analysis
We decided to quantify the conformational divergence of AMG 709 upon binding to the VEGFR-2 active site. For this purpose, aqueous biological medium was modeled in our ligand optimization procedure. Estimated binding energies for compound AMG 709 may be a direct outcome of varied internal energies of ligand in its protein bound and free states within biological media (ΔE instability ). ΔE inst. can  may account for the participation of solvation in binding profile. In the light of the above information, solvation energy of Motesanib molecule needs must be taken into account for the correlation of ΔE tb and ΔG b terms. This result might further demonstrate the important role of solvent molecules in determining final free binding energy of ligand-receptor system.
The estimated conformational change of ligand structure upon binding to the receptor was evaluated in a more detailed way via performing comparative conformational analysis of the molecular geometries. For this purpose, optimized 3D structure of AMG 709 was obtained by DFT calculations via B3LYP method in association with split valence basis set using polarization functions (Def2-SVP). Frequency calculation with same basis set was performed to confirm the optimized structure. All frequencies were real and no imaginary frequency was seen. The resulted geometric poses in terms of bond lengths and dihedral angles are summarized in Tables 2 and 3. It should be noticed that due to the uncertainty in the delicate position of hydrogen atoms in crystallographic file, associated data have not been shown in Tables. We found that all the calculated bond lengths of the DFT optimized structure were in adaptable correlation with the crystallographic data.
The varied dihedral angles between optimized and crystallographic ligand poses would be expected upon binding to the receptor active site. AMG 709 adapted some torsional distortions to get proper oriented pharmacophoric points. These well-oriented functional groups might be critical in achieving optimum key interactions with the residues of the VEGFR-2 active site.
Regarding the data in Table 3, some relatively significant angular deviations may be noticed. The observed rotation of C15-C16 bond ( Figure 5) let to the noticeable change in C8(13)-C15-N16-C18 dihedral angel (Table 3). This conformational distortion occurred at the amide linker.
All the mentioned conformational changes occurred in the structural moieties participated in interactions with key amino acids of VEGFR-2 active site (Figure 1). be defined as an energy difference for ligand in its free and protein bound states within aqueous medium. ΔE inst. needs to be considered in order to adjust obtained binding energies. Water was selected as a biological medium for this purpose.
For the purpose of calculating ΔE inst. , optimum structural conformation of compound AMG 709 was obtained in water and relevant energy was assigned to the free state . In the next step, the energy of receptor bound ligand was obtained in the crystallographic state. ΔE inst. may be well related with the free energy of binding via following equations: Higher ΔE inst. values support more positive total binding energies (ΔE tb ) consequently leading to weaker ligand-receptor interactions in terms of free binding energies (ΔG b ). Our calculations showed that AMG 709 tolerated 8.91 kcal/mol instability to gain the appropriate conformation in binding to the receptor. Based on the obtained results, ΔE tb was found to be -40.36 kcal/mol. Two conformational poses of the ligand are depicted in Figure 5.
However the difference between ΔE tb and ΔG b values associated with relevant ligand

Conclusion
Amino acid decomposition analysis provided further insight into the effect of individual amino acid residues on Motesanib/VEGFR-2 binding profile. Such structure-based studies may serve as efficient analyzing tools in evaluating the pharmacophore models. Owing to the prominent role of electrostatic forces in initial ligandreceptor interactions, charge-assisted H-bonds and cation-pi interactions need to be particularly attended. In the case of Cys919, the estimated binding energy could further confirm the role of this key amino acid in contribution to H-bond interactions reported for various VEGFR inhibitors. We could further demonstrate that Motesanib does not necessarily bind to the receptor in its optimum conformation state.